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Abstract 

The recent advances in alternating direct implicit (ADI) methods promise important new 

capability for time domain plasma simulations, namely the elimination of numerical stability 
limits on the time step. But the utility of these methods in simulations with charge and current 
sources, such as in electromagnetic particle-in-cell (EMPIC) computations, has been uncertain, 
as the methods introduced so far do not have the property of divergence preservation. This 
property is related to charge conservation and self-consistency, and is critical for accurate and 
robust EMPIC simulation. This paper contains a complete study of these ADI methods in the 
presence of charge and current sources. It is shown that there are four significantly distinct cases, 
with four more related by duality. Of those, only one preserves divergence and, thus, is 
guaranteed to be stable in the presence of moving charged particles. Computational verification 
of this property is accomplished by implementation in existing 3D-EMPIC simulation software. 
Of the other three cases, two are verified unstable, as expected, and one remains stable, despite 
the lack of divergence preservation. This other stable algorithm is shown to be related to the 
divergence preserving case by a similarity transformation, effectively providing the complement 
of the divergence preserving field in the finite-difference energy quantity. 
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1. Introduction 

Electromagnetic simulations of plasma can be made more robust with the use of implicit solvers, since 
these permit use of time step consistent with the phenomena of interest, rather than the fastest time-scale 
of the plasma, e.g., the plasma frequency, or speed-of-light Courant limit. An interesting class of those 
solvers is that of so-called Alternating Direction Implicit (ADI) solvers, of which one example was 
introduced in [1] and further analysis was presented in [2]. These solvers are found by separating the 
time development operator (curl acting on E and B) for Maxwell's equations into two parts. In each part, 
the derivative in any one direction (say x) acts on only one of the pairs (say Ey and B^) of the combined 
six components of the electric and magnetic fields. This allows one to develop a full update through 
Strang [3] splitting, which, upon repetitive application, alternates the updates of each of the operator 
parts separately. Furthermore, the simplicity of each part allows fully implicit updates with the 
requirement of only one-dimensional, tridiagonal solves, which can be rapidly effected through the 
Thomas algorithm [4]. 

In this paper we investigate how these methods might be combined with self-consistent charges and 
current, with the particular example being electromagnetic particle-in-cell (EMPIC) methods [5,6]. For 
such situations, the importance of charge conservation has been noted. A charge conserving current 
deposition algorithm has the property that integration of the finite difference continuity equation, 
A/7 = -V^^ • jAt (p is the charge density, j is the current density, and Vfd is the finite-difference form 
of the divergence operator), gives a charge density that is consistent with the charge density of the 
locations of the particles in the simulations. Such a current deposition algorithm was developed by 
Villasenor and Buneman [7] and later in [8] for smoother particles. As demonstrated in [9], failure to 
have a charge conserving current deposition algorithm can lead to nonphysical divergence buildup in the 
electric field that ultimately leads to catastrophic failure of the simulation. 

The second issue, which is dealt with in this paper, is whether the electromagnetic update algorithm 
preserves this divergence of the current. Should indeed this be the case, then one has the highly 
desirable property that, if V^^, - 8 = and V^^ -E = yo/^g at one time-step, then these remain true, to 
machine precision, after application of the update operator, to the next time step. Stated another way, 
for the magnetic update, the requirement is that a divergence-less field, upon update, leads to a 
divergence-less field. For the electric update, the associated condition is that, given a conservative 
charge deposition scheme. Gauss's Law, V^^ AE = -V^^ • )At/ , remains valid, to machine precision, 
after application of the update operator, without the need for any additional divergence cleaning step. 
Such divergence preservation is indeed the case for the Yee [10] update and for a Crank-Nicholson 
update as well, but it is not a priori clear for the ADI updates. 

In this paper, we perform the splitting of the curl operator in Maxwell's equations into two parts, 
denoted M and P, leading to the four fundamental operators of ADI algorithms, which we denote as 
(l + f m), (l + f p), (l-f m)"^ and (l-f P)"\ The first two operators are explicit, and the second 
two are implicit, since they are inverses. An ADI algorithm uses all four of these operators to construct 
the full update sequence. There are 24 = 4! possible combinations of the four operators. Of these 24, 12 
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are equivalent to each other by simple duality, that is, interchange of M and P. Of the remaining 12 
combinations, 4 are not transformable to Strang splitting, because successive alternation between M and 
P operators prevents construction of unitary operators. Finally, of the remaining 8 combinations, only 4 
are time-reversible, which is an important property of the Maxwell update and assures a minimum of 
2"''-order accuracy. Of those four combinations, one arises naturally from that introduced by [1] but 
does not preserve divergence. Of the other three combinations, only one is found to be divergence 
preserving. Strikingly, this one divergence preserving ADI composition has not been previously 
identified or studied. Finally, the remaining two compositions are related to the previous two by 
similarity transformation, and as previously indicated, are not divergence preserving. 

We then proceed to examine the consequences of using any of these updates in EMPIC, having 
implemented them in the VORPAL [11] computational application. We show that, like the previous 
results of [9] for charge non-conservation, the lack of divergence preservation in two of the algorithms, 
including the one most closely related to that of [1], leads to catastrophic failure of the associated ADI- 
EMPIC simulations. We show that the new, divergence preserving algorithm does not have this 
catastrophic failure. Nor does the algorithm related to it by a similarity transformation. 

These results can be applied to ADI algorithms in simulations of many phenomena for which one would 
like to use EMPIC, but for which the Courant limit is constraining. One example comes fi-om the 
requirement to resolve the plasma Debye length in order to avoid self-heating [12]. This requirement in 
an explicit EM simulation then leads to a Courant time step that is A? « /l^ / c . However, if light waves 
are not important, one really needs a time step that can be much larger, of order the time for a thermal 
electron to cross a cell, or A/^ w /l^ / , where v^, is the electron thermal speed. Thus, use of explicit 
methods, as might be needed to capture magnetic effects, requires a much smaller time step and, hence, 
much greater computational effort. Another common situation involves the propagation of a bunched, or 
otherwise spatially varying charge profile beam through a tube or cavity. The beam ultimately comes to 
equilibrium, with its self electric and magnetic forces in balance with any external fields and the beam 
divergence. Now the time scale is very long, as one is computing an equilibrium, so the Courant 
stability condition is even more restrictive. In general, the Courant condition is constraining whenever 
one is dealing with systems where one must keep magnetic effects yet light waves are not important to 
the dynamics. 

The organization of this paper is as follows. In the following section we review the ADI-EM algorithms 
and derive the four fundamental operators fi-om which all possible update operators can be assembled. 
We then show, in Sec. 3, using time-reversal arguments, that there are two possible second-order update 
operators for vacuum electromagnetics, and even those are related by a similarity transformation. In Sec. 

4, we introduce current sources. This breaks the equivalence, so that there are four update operators, 
consisting of two pairs, with each in a pair related to the other by a similarity transformation. We show 
that only one of the four operators is divergence preserving. In Sec. 5 we present numerical results for 
ADI-EMPIC. These results show that only the divergence preserving operator and its relative by 
similarity are stable; the other two algorithms lead to catastrophic failure in ADI-EMPIC simulations. In 
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Sec. 6 we derive the steady state solution for these two operators and show that the one similar to the 
divergence preserving update operator has the usual finite-differenced steady state Maxwell equations. 
In Sec. 7 we derive an energy invariant for the divergence preserving operator and its relative by 
similarity, and we show how its energy differs from the standard divergence preservation of FDTD EM. 
In Sec. 8 we discuss the use of these algorithms in conjunction with complex boundaries. The last 
section contains a summary and some conclusions. 

2. ADI EM algorithms 

The ADI algorithms derive fi-om the fact that Maxwell's equations, 
dB 



dt 



= -VxE 



(la) 



and 



8E 

~dt 



+ c^VxB, 



(lb) 



can be written in the form. 



dv 
dt 



=s+^+m)v 



(2) 



where V is the six-component field, (E, cB), S is the six-component field (-j/so, 0), and the operators 
P and M are defined by 



P V = P 







dcB^ jdy 






dcBjdz 






dcBy 1 dx 




= c 








dE^ 1 dx 


SB._ 







and M-V = M' 







~dcBy/dz~ 


Ey 




dcBjdx 


E. 




dcBjdy 


cB, 


= -c 


dEjdy 


^By 




dEjdz 


yB._ 




dEy jdx 



(3) 



The mnemonic is that the operator P has a plus sign on the right, while the operator M has the minus 
sign. As can be seen, the interchange of P and M operators is equivalent to an interchange of E and cB, 
together with sign (parity) exchange. In the remainder of this paper we will assume what might be 
called an "M-first" choice of ADI-duality, and simply state here that the analysis proceeds identically 
under the interchange of P and M operators. 

For numerical differencing, the Yee layout of the fields is assumed. In this layout, the electric fields are 
centered at the edges of a cell (indexed by ij,k), while the magnetic fields are centered at the faces of the 
cell. We adopt the following notation to indicate the transition from field representation to discrete 
spatial representation: the previous tilde quantities are fields and calculus-operators, whereas without 
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the tildes, the quantities are arrays and finite-differencing matrices. Thus, the finite difference forms of 
the above operators are 

c{b ..-B . . , , )/Az" 
(e . -E . . J/Az 

x,ij ,k+\ x,i,j,k // 

(e -t-E ■ ■.)/Ax 



P V = c 

















-By,- 























and M • V = -c 



(4) 



where now V is the full array of all values for all components and all cells, and P and M are linear 
operators on that space. Thus, the finite differenced Maxwell equations, in the absence of sources, are 



5V 

dt 



= (p + m)-v 



(5) 



An important property of the discretized Maxwell system of equations is that the matrix operators are 

T T 

anti-symmetric, that is, M=-M and P=-P , where the superscript T' denotes transpose. 

With Strang splitting, we break the above equation into two, and we solve each one separately, giving 
two advance operators. If each operator is found to second-order accuracy (third-order error), the full 
update can be found by application of the square root of the first operator followed by the second 
operator and then the square root of the first operator. More definitively, we first consider the equation. 



8X 



= P X 



dt 

which becomes 



X"^' -X" =f P-^"^' +X"] 



(6) 



(7) 



upon 2"''-order time-centered finite differencing, at discrete time-steps, t„=nAt using the usual time- 
superscript notation, X"=X(t„). The update solution of this equation is, of course, 

X"-"' =Tp -X" =(l-f P)"' -(l + f P)-X" , (8) 

thus defining the unitary second-order accurate time-advance matrix operator, Tp, based upon the 
original P operator. Similarly, 



T^=(l-fM)-'.(l + fM) 



(9) 



gives the second-order accurate time-advance based upon the M operator. The unitary property of these 

T 

matrices-operators, that is, Tp Tp=l and similarly for Tm, is a result of the anti-symmetric property of 
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the original P and M matrix-operators, and the commutability of the matrix factors, ^ + ^m)^' and 
(l - ^m), and analogously for P. Consequently, following Strang splitting, the composite operator, 

Strang " ' V^"/ 

is also second-order accurate. 

This is not the update proposed by [1]. Instead, they proposed the ADI-update operator, 

Ti =(l-f M)"'-(l + f P)-(l-f P)"'-(l + f m), (11) 

which, we show in the next section, is also second-order accurate, although this property was not fully 
advertised in [1]. The two operators are related by a similarity transformation. 

T,=(l-^M^p.T,,,„^.(l-4M^r. (12) 

Thus, the spectrums of the two updates are identical, and the update of one corresponds to the update of 
the similarity transform applied to the state vector of the other. Lee and Fomberg [2] noted that the 
operator of [1] could be rewritten, 

T^,,=(i-fMy.(i-fpy.(i+fp>(i+fM), (13) 

by virtue of the commutability of the middle two operators in Eq. (11). Application of this operator in 
an update equation for V allows the inverses to be taken in turn, producing a traditional ADI algorithm, 

y«+l_y„ ^^(M + p).(y«+l+y«^_Mp.M.(y«+l_y«^ (14) 

that differs from the second-order Crank-Nicholson operator by the last term that is 0(At ). Hence, the 
update operator (13) is second-order accurate. The notation adopted for this operator, "GSS", refers to 
"Gradient Steady State," due to the fact that this update, Eq. (14), allows a steady-state condition, 
y «+i ^ y« ^ whenever V" is a pure gradient, since the gradient is in the null-space of the curl operator, 
(m + p), for Yee-cell finite-differencing. This important property will be explored in more detail in 
Section 6. 



3. Second-order accurate ADI operators 

Inspired by the results of Lee and Fomberg, we approach this from another direction. Namely, we 
consider any and all possible operators that can be constructed from a product of each of the four 
operators appearing in Eq. (11), but we restrict ourselves to operators that are time -reversible. Time 
reversibility guarantees second-order accuracy, because time reversibility guarantees that the first non- 
vanishing term in the power series is third order in the time difference, just like in the actual evolution, 
and so the first possible difference is O(At^). 
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First we must understand how time reversal acts on the operators of Eq. (11). Since time reversal is 
computing the final state in terms of the initial, the order of the operators must be reversed, and the 
inverses interchange with multiplication. For an operator to be time reversal symmetric, then its inverse 
followed by /S.t —>-At must give the same operator. We see that the operators (11) and (13) introduced 
in Refs. 1 both have time-reversal symmetry and are, therefore, second-order accurate. 

We now enumerate the other operators that can be time-reversal symmetric. As noted previously, we 

restrict ourselves to one specific choice of the P M interchange duality, namely that all sequences of 
the form of Equation (11) can be taken to begin with either (I-^mJt' or (i + ^m) and end with the 
other corresponding time reversed operator. Thus, we have the previously defined update-operator from 
[1], and the operator that results from exchanging its first and last terms, 

Ti =(l-fMy -(l + f P)(l-f Py -(l + f m), (15a) 

and 

T2 =(l + f M)(l + f P)(l-f P)r' -(l-f My . (15b) 

Each of these two update operators has a corresponding equivalent, in the absence of currents, with the 
center terms commuted, including the previously noted Eq. (13). 

T^,,=(l-fMy.(l-fPy.(l + fP>(l + fM) (15c) 



T^^ = (1 + f m) (1 - f Py • (1 + f P) (1 - f My ( 1 5d) 

The notation adopted for this last operator, "Z)/'", refers to "Divergence Preserving." This previously 
unheralded update operator will be the subject of the next section and is indeed the impetus for this 
study. 

These two pairs of operators are related by similarity transformation. 



T 



2 



^1-^M^ |-Ti 



V^M^V and Tj^pJl-^mA-Tcss-il-^MA'. (16) 



4 



V 4 



V 4 



V 4 



Thus, in the absence of currents, there is only one fundamental second-order accurate vacuum 
elecfromagnetics update operator, with all others related to it by duality, commutability, or similarity 
transformation. Because the operator appearing in Eq. (16) arises repeatedly, we denote it as R, and 
discuss it further. Note that, in the continuous field representation, is a component- wise diagonal 
operator, which can be represented in derivative and specfral form as. 
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d'cBjdy' 
d^cBjdz' 
d^cBjdx' 



klcB^ 
klcB^ 
klcB, 



(17) 



so that 















































where for continuum fields, the spectral-form coefficient is 



kfibdil 









l+a! 



(18) 



(19) 



For discrete representation, the matrices, M and R are the block tridiagonal matrices based upon the 



2"'^-deriviative finite-differences, and the spectral-form coefficient is 



ctst . 

(J^ = Sin 

Ax 



A:,. Ax,. 



V 



(20) 



The R matrix is positive definite. For well resolved variations, kiAxi,«n and kicAt«2, it is 
approximately unity to second order, and departs from unity for poorly resolved variations in proportion 
to the Courant ratio, cAt/Ax. Finally, we note that R contains only the M operator and not the P operator. 
Thus, R serves as an indicator of the choice of the ADI duality. 



4. Divergence preservation 

Divergence preservation becomes an issue when current sources are added into Maxwell's equations. In 
EMPIC, the electric current, j, is usually computed at times halfway between the discrete times at which 
the electric field is known. Consequently, charge, p, is computed at the same times as the electric field, 
enabling Gauss's Law, with its divergence operation, to be evaluated. For ADI, we would like to 
continue to add the current in a time-centered manner, and in the simplest way possible. Thus, let us 
assume a current source evaluated at the half time step, S""^*^". Thus, we look at possible update 
algorithms where the current is added once, in the center of the operators (15), e.g.. 
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v;^' = Ui (v; ) = (i - f m)"' • (i + f p)- [(i - f p)"' • (i + f m)- Vi" + A^s"^^ _ 

Vf' = U2 (V^" ) = (1 + f m)- (1 + f P)- [(1 - f P)"' • (1 - f m)"' • \l + MS"^^' _ 

=u^,,fc)-(i-f m)" .(i-f p)-' .[(i+f p).(i+f m).v^,, +A/S"^>^], 



and 



V;;^=U^pfc)-(l + fM).(l-fp)-^[(l + fp).(l-fM)-'.V^^+A/S 



(21a) 
(21b) 
(21c) 

(21d) 



With this introduction of a current source term, the updates that were equivalent by virtue of 
commutation are no longer equivalent, since the source is added between the two commuting terms. 
Hence, there are four distinct second-order updates when current is present, rather than two, plus four 
more that can be obtained by duality. 

In the Yee algorithm, divergence preservation comes from the fact that, as in the vector identity, 
V-VxA=0, the numerical curl operator, (P+M), is in the null-space of the numerical divergence, which 
we denote here by V^^ • , and note that this provides the separate curls on the E and cB parts of the field. 
Thus, for any field array, X, 



V^^-(P + M)-X = 0. 



(22) 



We now analyze the divergence preservation properties of Eq. (2 Id), \Jdp, on a field \dp, for which 
purpose we rewrite as 



(i-f p)(i+f M)r' -v;;;' =(i+f p)(i-f m)"' ■y"^p+Ats" 

Divergence preservation follows from the following identities, 

(1 - f p) (1 + f m)"' = (1 - f p - f m + f m) (1 + f m)"' 

= l-f (P + M)-(l + f m)"' 

and, similarly, 

(i+f p)(i-f M)r' =i+f (p+M)(i-f M)r\ 



(23) 



(24) 



(25) 



Plugging Eqs. (24-25) into Eq. (23), and collecting terms results in an analogue to Eq. (14), although it 
is not as immediately recognizable as an ADI update. 



Kp -v;;^ =f (M+p)-[(i+f m)"' -v;;;' +(i-f m)"' •v;;^J+a^s 



(26) 



We take the divergence of both sides of Eq. (26), and note from Equation (22) that the divergence, 
V^^ • , operating on (M+P) vanishes, thus yielding. 
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(27) 



which is precisely the divergence preservation property, including source current. Since there are no 
magnetic currents, this result guarantees that the numerical divergence of the magnetic field always 
vanishes if it did so initially. Similarly, it implies that, provided the current and charge density satisfy 
the numerical continuity equation, the numerical divergence of the electric field will always equal the 
charge density if it did so originally. That is, 

V^^E -E;;' -V^^, ■El,=-V,^^-Atr"'/^o=(p"'' -p"ho- (28) 

where Vfde- is the part of the divergence matrix associated with the electric field. The above 
manipulations relied on the fact that the first applied operator was an inverse, e.g., implicit, while the 
second was explicit, so that one has the identity plus a term whose left-most term is precisely the curl 
operator. Remarkably, of all the updates in (21), only Vdp, (2 Id) has this form. Hence, none of the 
other updates will preserve divergence (though the dual obtained by interchanging P and M will have 
the same divergence preservation). All four updates have been implemented in numerical software to 
confirm this unique divergence preserving property, and also to verify the lack of divergence 
preservation in the other three updates. This will be discussed in fiirther detail in Section 5. 

In the absence of current sources, similarity transformations by the matrix R related these update 
operators, l<-^2 and GSS<r^DP. Thus, we expect that these pairs will have similar stability properties in 
the presence of currents and charges. In particular, while the electric and magnetic fields for the update 
Vgss are not updated in a divergence preserving way, they are related by a transformation to electric and 
magnetic fields that are updated in a divergence preserving way, e.g., to the electric and magnetic fields 
for the Udp update. More precisely, the DP and GSS solutions are related by 

V^P=R-V^,,. (29) 

Now, since the \dp solution is divergence preserving, it implies that the GSS update preserves not 
divergence of the field, but rather divergence of the field multiplied by the R matrix, e.g., denoting 
VpDE- and Re to be those parts of the divergence and R matrices associated with the electric field, we 
have, 

V^z,E • Re • E^;^ - V^^ • R, • E^,, = -V^^ • Atr'" I = (p"^' - p-y^, , (30) 

This operator, V^^j, -Rj. •, has a larger stencil than the V^^j. •. Also, as the spectral representation of 
the R matrix, Eq. (18), enhances larger wavenumbers, it acts as an anti-smoothing operator. In addition, 
as this operation, V^^j, • - ^css ■> resembles the use of non-trivial dielectric, 8, in Gauss's Law, V-e-E, 
it is natural to think of the R matrix as an additional numerical permittivity/permeability associated with 
the ADI finite-differencing algorithm, especially if one recalls that R departs from the identity in 
proportion to worsening resolution of the field variations, and that R is an indicator of the choice of ADI 
duality. Further exploration of the role of the R matrix is provided in Sec. 6 and 7. 
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5. Numerical results for self-consistent charges and currents 

To determine the effect of divergence preservation we developed prototype implementations within the 
VORPAL [11] computational framework, which can be used for EMPIC. We chose parameters to 
match the previous [9] numerical results. Namely, our simulations were 2D inside a box, Im on a side 
and with 20x20 cells. The simulation was initialized to have no particles or fields, but then 5A of 30 
keV electrons were injected from the middle third of the left wall. The time step was chosen to be the 
Courant value, and all four algorithms in (21) were tested. The simulations were run for 80,000 time 
steps, which is slightly over 900 times of transit of the beam across the system. 

Figure 1 shows the results for the x-y scatter plot of the beam after 725 transit times using the ADI 
update operator, Uj , based upon [1], with currents added in the center. One can see that the beam is 
developing an unphysical divergence, much like that seen in [9], when a non-charge conserving current 
deposition algorithm was used. While this is for one instant in time, further integration shows that the 
situation gets worse. It is accompanied by a growth in the divergence error, V^^ •£-/?/£■(,, which 

6 2 

vanishes initially, but by this time in the simulation has grown to peak values of 3x10 V/m . By 
comparison, the charge density for the beam is about 7.5x10"^ C/m^, or pi s^^ « 8.5e3 V/m^. Thus the 
error in the divergence, by this time, is orders of magnitude greater than the actual beam charge density. 

In contrast. Fig. 2 shows the results for the x-y scatter plot of the beam after 725 transit times using the 
newly introduced divergence preserving ADI update operator, \J£,p . There is no sign of beam 
divergence developing, as expected. The behavior is reminiscent of the standard Yee update with 
charge conserving current deposition that was shown in [9]. For this case, the divergence error is around 
3x10 V/m , or around round-off error when compared with pi Sq ^ 8.5e3 V/m . 

As noted earlier, the field solutions from update operators Uj and U2 are related by simple 
transformation. Hence they are expected to have identical stability properties in the presence of particles. 
Indeed, simulations show this to be the case. Use of the operator also leads to unphysical beam 
divergence after a few hundred transit times. 

Similarly, the field solution using U^^^ is stable, as it is related to field solution of V^phy a simple 
transformation, and simulations using U^^^ are seen to have the same basic properties as those using 
\}p,p, i.e., no nonphysical beam divergence. For these simulations using U^^^ , we observe a divergence 
error that is on the order of 300 V/m , which is less than, but not negligible, compared with the beam 
value of pi ^ 8.5e3 V/m . So even though the usual divergence is not preserved by the operator, 
Vqss ' f^'^t th^t ^ related divergence quantity is preserved appears to be enough to assure stability. 

6. Steady state solutions 

Because the update U^i^ preserves the divergence, we are guaranteed that a solution using this operator 
always has V^^ ■ B = and V^^, • E = /j/^-q , in exact analogy with the continuous-field representation 
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of Maxwell's equations. However, the traditional Yee-cell finite-difference electromagnetics has the 
further property of steady-state solutions satisfjdng V^^xE = and V^^xB = //oj, also in analogy 
with the Maxwell's equations. Here we investigate the steady state solutions for some of the previously 
defined ADI updates. 

We are particularly interested in the steady-state solution of the divergence preserving update, Eq. (26). 
Setting V^p , V^p V^p in this equation, and noting that 

(l + f m)-' -V;;' +(l-f m)-' .v^p ^[(l + f m)-^ +(l-f m)-^} V^p =2R-^ • V^p (31) 
we can write the steady-state solution of the divergence preserving update in the form, 

-(P + M>R-' -V^p =S, (32a) 

Noting the previously established relationship between Ydp and Ygss, Equation (29), it follows 
immediately that 

-(P + M>V^,,=S, (32b) 

which is the analogue to the continuous-field Maxwell equations in steady-state. Thus, the fields found 
from using U^ip do not satisfy the normal steady-state equations, rather they are related by the R matrix 
to the fields, Vgss, that do. Because the R matrix enhances shorter wavelengths, this ADI solution will 
have shorter wavelengths enhanced. To summarize, the GSS solution has analogous Maxwell steady 
state behavior, but it has divergence error. In contrast, the DP solution has no divergence error, but it 
does not satisfy the usual finite difference Maxwell equations in steady state. 

It should be noted that this study has investigated the ADI algorithms for Maxwell equations in vacuum. 
After demonstrating the properties of the DP and GSS solutions. Equations (28) and (32b), we would be 
remiss not to point out the obvious similarity between the vacuum ADI algorithm and the Maxwell 
equations in the presence of non-trivial materials. When materials are present, the non-tiivial dielectiic 
creates a role separation between electric field, E, and electric displacement, D=8-E, and similarly with 
B=(j,-H. In this role separation, it is the fields D and B which occur in the divergence equations, and it is 
the fields E and H which are acted upon directly by the curl operator, and thus occur in the steady-state 
equations. The DP and GSS vacuum solutions have similar role separation, namely \dp occurs in the 
divergence equation, and Ygss occurs in the steady-state equations, and furthermore, as already 
mentioned, the R matrix connects Vdp to Vgss exactly as e connects the D to E and \i connects B to H. 
Thus, despite the fact that we are talking specifically of vacuum field solutions, it is not a difficult 
observation to note that the DP solution behaves analogously to D and B in a material, while the GSS 
solution behaves analogously to E and H in a material. Whether this observation is a useful artifice 
remains to be seen. More to the point, we leave the investigation of the ADI algorithms in the presence 
of actual materials for later study. 
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7. Energy conservation 

We consider only the GSS and DP updates, as those are the only cases that are stable in the presence of 
charged particles. In [2] it is shown that Toss has a positive energy quantity, which corresponds to 



W = V 



GSS 



R- V = V 

^ GSS ^ DP 



R 



■ V = V 

^ DP ^ GSS 



DP ' 



(33) 



where the 2"^ and 3'^'' definitions follow from Eq. (29). This quantity is conserved, as is the case in that 
analysis, in the absence of charges and currents. The goal of this section is to determine how this 
quantity evolves in the presence of charges and currents. We seek the change in the vacuum energy 
quantity (33), which can be written in the form. 



AW 



Vn+l -« 7/7+1 



-V" -V" 

' GSS ^ DP • 



(34) 



This quantity arises naturally if one multiplies each side of Eq. (26) by the term following the curl, 
noting that due to the anti-symmetric property of the curl, that term vanishes, hence leaving just 



A?s""^^ • [(i + f m)"' • v;;' + (i - f m)"' • v^p] 



(35) 



Substituting from Eq. (26) for the left values of V^p gives. 



- V^,, ) • [(1 - f M ) • V - ' + (1 + f M ) • V^, J = 

Ats"^^' ■ [(i + f m)"' • v;;;' + (i - f m)~' • v^^] 



(36) 



Finally, we note that by virtue of the anti-sjTiimetric property of the matrix M-R, the quantity 
fe - nss ) M • fe^ - V«p ) is zero, leaving just. 



H+l 



AW = V • V 

^ GSS ^ DP 



■ V" -V" - AtS' 

^ GSS ^ DP ■"'^ 



(l + ^M)-'.V^;'+(l-fM)-^V^,J 
= A^S"^>^ -[(l-f m). V«J^ +(l + f m). V^,,] 



(37) 



Equation (37) is similar to the expression for mechanical energy in the usual explicit finite-difference 
Maxwell equations, with the novel feature being the presence of the factors that make up the R matrix. 
We see from Eq. (33) that since the R matrix acts as an anti-smoothing operator, the GSS field must 
contain less short-wavelength field amplitude than the DP field. Similarly, the ADI form of the 
mechanical energy is smoothed if expressed in terms of the DP fields, and anti-smoothed if expressed in 
terms of the GSS fields. A full study of fluctuations in this system is a subject for future research. Here 
we note that if there is a fixed amount of energy in each Fourier mode, as one might expect from 
equipartition, the electric and magnetic fluctuations will be larger for the DP case, as the mode energy 
has a smaller factor in front of it. 
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8. Boundary conditions 

Usage of these operators with boundary conditions is straightforward. Typically, two different types of 
boundary conditions are used. Stair-step boundary conditions are those for which a whole cell is either 
included or excluded from the simulation region. Cut-cell boundary conditions [13,14], are those for 
which the magnetic update is modified through use of Faraday's law for the fraction of the face that is 
within the simulation region. 

For stair-step boundary conditions, with a cell either entirely included or excluded from the simulation 
region, the boundary ultimately consists of the faces of all included cells. Boundary conditions are 
applied by setting the values of the electric field on all cell edges on the boundary. The magnetic update 
equations are then valid for all interior faces as well as those on the boundary. The matrices P and M 
are then found by setting to zero any matrix elements operating on an electric field for which the 
corresponding edge is exterior - e.g., not even on the boundary. These terms give the time derivatives 
of the magnetic field. For the complement, one can use the regular magnetic update everywhere, 
provided it is followed by setting the values of the electric field on all boundary edges to the boundary 
values. 

For the Dey-Mittra type cut-cell boundary conditions [13], which provide second-order accuracy in 
global quantities, again the electric update is unchanged, except that the electric field for any edge 
wholly outside the simulation region is not updated but remains at the boundary value. However, the 
magnetic updates are modified by changing the coefficient of the electric fields in P and M by 
multiplying them by the relative length of the edge within the simulation region divided by the area of 
the face of the magnetic field that is within the simulation region. For faces entirely outside the 
simulation region, the corresponding elements of the P and M matrices can be set to zero. 

The Dey-Mittra boundary conditions are known to lead to a reduction in the maximum stable time step, 
when traditional leap-frog updating is used. The uniformly stable update algorithm [14] eliminates this 
reduction. We note here the important fact that use of the ADI methods, in lieu of leap-frog, also 
eliminates this time-step reduction, since the ADI methods are stable for arbifrary time step. 

9. Summary and conclusions 

An exhaustive study of ADI update operators for elecfromagnetics in the presence of charges and 
currents has been completed, and all possible 2"'*-order ADI update operators having time-centered 
current addition at only one place have been identified. Of these only our newly introduced divergence- 
preserving U^p update from (2 Id) and its dual (from interchanging P and M) are divergence preserving. 
We also demonstrated the curl-steady-state property of the ADI algorithm of [2], generalizing it to 
include current, and resulting in the l^ass update from (21c), and demonstrate that a simple 
transformation connects the divergence-preserving and curl-steady-state field solutions. Upon 
implementation in EMPIC software, it was observed that only U^.^ and U^^^ are suitable for use in 
particle simulations, as the others show long time divergence error growth that leads to unphysical 
behavior. 
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Figure 1 
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Legend: Configuration-space scatter plot of a beam after 575 transit times for the update operator Ui. 
Artificial charge build-up on the grid causes an unphysical beam divergence instability. 
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Figure 2 
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Legend: Beam after 775 transit times using the charge conserving update, \Jdp- Divergence error is 
zero to machine precision, and so beam transport is stable. 
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Figure Legends 



Figure 1. Configuration-space scatter plot of a beam after 575 transit times for the update operator Ui. 
Artificial charge build-up on the grid causes an unphysical beam divergence instability. 

Figure 2. Beam after 775 transit times using the charge conserving update, Vdp- Divergence error is 
zero to machine precision, and so beam transport is stable. 
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